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The adiabatic insertion of a n flux into a quantum spin Hall insulator gives rise to localized spin 
and charge fluxon states. We demonstrate that tt fluxes can be used in exact quantum Monte Carlo 
simulations to identify a correlated Z2 topological insulator using the example of the Kane-Mele- 
Hubbard model. In the presence of repulsive interactions, a tt flux gives rise to a Kramers doublet of 
spin fluxon states with a Curie law signature in the magnetic susceptibility. Electronic correlations 
also provide a bosonic mode of magnetic excitons with tunable energy that act as exchange particles 
and mediate a dynamical interaction of adjustable range and strength between spin fluxons. tt 
fluxes can therefore be used to build models of interacting spins. This idea is applied to a three-spin 
ring and to one-dimensional spin chains. Because of the freedom to create almost arbitrary spin 
lattices, correlated topological insulators with tt fluxes represent a novel kind of quantum simulator 
potentially useful for numerical simulations and experiments. 
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I. INTRODUCTION 

A topological insulator represents a novel state of mat- 
ter characterized by a special band structure that can re- 
sult, e.g., from strong spin-orbit interaction [UH]. In two 
dimensions, this state is called a quantum spin Hall insu- 
lator, and has deep connections with the quantum Hall 
effect, including the coexistence of a bulk band gap and 
metallic edge states, the absence of symmetry breaking, 
and the possibility of a mathematical classification [31 [4J . 
Importantly, because of the absence of a magnetic field, 
the quantum spin Hall insulator preserves time-reversal 
symmetry which provides protection against interactions 
and disorder [SHZ] . The quantum spin Hall insulator has 
been realized in HgTe quantum wells [HI [5] . 

Correlated topological insulators with strong electron- 
electron interactions are the focus of current research |10j . 
Intriguing concepts include electron fractionalization in 
the presence of time-reversal symmetry [TTHT3] , spin liq- 
uids P3HH], an d topological Mott insulators [T7JQ2]. Re- 
markably, some of the theoretical models can be studied 
using exact numerical methods. A central problem in 
this context is the question of how to detect a topological 
state directly from bulk properties, for example in cases 
where the bulk-boundary correspondence breaks down. 
Experimentally, this issue also arises in the absence of 
sharp edges in proposed cold atom realizations as a result 
of the trapping potential [19] [20] . The classification in 
terms of a Z2 Chern-Simons index relies on Bloch wave- 
functions and is therefore only valid for noninteracting 
systems. Generalizations involve twisted boundary con- 
ditions [21] or Green functions [22~f-i27|. and are challeng- 
ing to use in experiments or exact simulations. Indirect 
signatures such as the closing of gaps [TB] or the cross- 
ing of energy levels |28j require, among other difficulties, 
experimental tuning of microscopic parameters. 

Topological insulators show a unique response to topo- 



logical defects such as dislocations [55J ISO] or tt fluxes 
[l2l 130] [31] . Upon adiabatic insertion of a tt flux, Fara- 
day's law together with the quantized transverse conduc- 
tivity gives rise to midgap charge and spin fluxon states 
[T2l [31] . These states are exponentially localized around 
the flux [121 [31]. The existence of these states is en- 
sured even in the presence of interactions or disorder 
by time-reversal symmetry, and has been suggested as 
a bulk probe of the Z2 index [T2J [31]. The concept of 
fluxons can also be generalized to situations where spin 
is not conserved, such as in the presence of Rashba cou- 
pling. In three dimensions, a magnetic flux gives rise to 
the wormhole effect |32j . Electron-electron repulsion lifts 
the degeneracy of charge and spin fluxons, but the two 
degenerate spin fluxon states constitute a localized spin 
with S z — ±1/2 [T2]. Dynamical tt fluxes emerge in the 
context of fractionalized topological insulators P21 US] . 

Previous work on tt fluxes in noninteracting quantum 
spin Hall insulators (TSJ [30] [31] was based on square- 
lattice models such as that for HgTe quantum wells [5J. 
Here we consider the half-filled Kane-Mele model on the 
honeycomb lattice [3], historically the first model with 
a Z2 topological phase, which has close connections to 
graphene [3J, the integer quantum Hall effect [33], and 
when including interactions, to correlated Dirac fcrmions 
[15j . Topological phases of interacting fcrmions on hon- 
eycomb lattices may be realized in transition metal oxides 
[33] , semiconductor structures [35] , graphene [33] , or cold 
atoms [37], see also Ref. [TU] , 

Here we use tt fluxes in combination with exact quan- 
tum Monte Carlo simulations, and show that they can be 
used efficiently to probe the topological invariant of cor- 
related topological insulators. In particular, this method 
does not rely on an adiabatic connection to a noninter- 
acting state, and may also be used for fractional states. 
In addition, we demonstrate that tt fluxes permit the 
construction of quantum spin models of almost arbitrary 
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FIG. 1. (Color online) For a lattice with periodic boundaries, 
7r fluxes can be inserted in pairs. Each flux threads a hexagon 
(highlighted in blue) of the honeycomb lattice, and the pair 
is connected by a branch cut (blue line). Hopping processes 
crossing the branch cut acquire a phase e 17T — — 1. 



geometry and with tunable, dynamical interactions. The 
spins correspond to the spin fluxons created by insert- 
ing 7r fluxes, and the interaction is mediated by mag- 
netic excitons corresponding to collective magnetic fluc- 
tuations of the topological insulator. These spin models 
can be studied theoretically with the quantum Monte 
Carlo method, or experimentally. As examples, we show 
that a ring of three spins has a ground state with magne- 
tization 1/2, and that a one-dimensional chain of fluxons 
undergoes a Mott transition and is described at low en- 
ergies by an XXZ model. 

The article is organized as follows. In Sec. [TTJ we in- 
troduce the Kane-Mele and Kanc-Mclc-Hubbard models. 



Section [TTT| provides details about the methods. The use 
of 7r fluxes as a probe for topological states is discussed in 
Sec. |IV[ whereas the construction of quantum spin mod- 
els is the topic of Sec.|V] Conclusions are given in Sec.|VI[ 
and we provide three appendices. 



II. MODEL 

The half-filled Kane-Mele model with additional 
electron-electron interactions can be studied with pow- 
erful quantum Monte Carlo methods [HI [35] . Using the 
spinor notation cj = (c^, c]^) , where c| CT is a creation op- 
erator for an electron in a Wannier state at site i with 
spin a, the Hamiltonian reads 

H K m = ~t^2 T v&o + 1 A T » ^ ("v ' cr ) S (*) 

<i.i> ({id)) 

+ \a y £ j T ij c\(s x d tj ) -zcj. 

(i,j) 

The notation and indicates that the sites i 

and j are nearest neighbors and next-nearest neighbors, 
respectively, and implicitly includes the Hermitian con- 
jugate terms. 

The first term describes the hopping of electrons be- 
tween neighboring lattice sites. The second term repre- 
sents the spin-orbit coupling which reduces the SU(2) 
spin rotation symmetry to a (7(1) symmetry. The third 



term is an additional Rashba spin-orbit coupling |39j . 
The additional factors — ±1 take into account any 
7r fluxes present, whereas the original Kane-Mele model 
(without 7T fluxes) is recovered from Eq. ([T]) by setting 
Uj = 1. 

The spin-orbit term corresponds to a next-nearest 
neighbor hopping with a complex amplitude iA, and has 
been derived from the spin-orbit coupling in graphene 
[3J. This hopping acquires a sign ±1 depending on its 
direction, the sublattice, and the electron spin. This sign 
is encoded in (i/^- • cr), where 



dik x d 



\d lk x d, 



kj\ 



(2) 



dik is the vector connecting sites i and fc, and k is the 
intermediate lattice site involved in the hopping process 
from i to j. For a coordinate independent representation, 
the vectors d a p are defined in three dimensional space, 
although the z component vanishes. The vector cr is 
denned by cr — (a x ,a y , a z ), with the Pauli matrices a a . 

The last term in Eq. ([!]) is the Rashba spin-orbit in- 
teraction [3J [3] . It is defined in terms of the spin vector 
s = er/2, and the unit vector dy which can be expressed 
in terms of the nearest-neighbor vectors 6\, 62, S3 [40] . 
The Rashba coupling breaks the z 1— ¥ —z inversion sym- 
metry, and has to be taken into account, for example, in 
the presence of a substrate. Because this term includes 
spin-flip terms, spin is no longer conserved. The Rashba 
term has been included in the results for the noninteract- 
ing model ([I]) , but cannot be included in quantum Monte 
Carlo simulations of the interacting model ^ due to a 
minus-sign problem. 

The model can be solved exactly [31 El ST]. In 
the absence of Rashba coupling, a — 0, the Kane-Mele 
model describes a Z 2 quantum spin Hall insulator for any 
A > 0. This state is characterized by a bulk band gap 

A sp = 3^3 A, a spin gap A s = 2A sp , and a quantized 

2 

spin Hall conductivity a s xy = The topological state 
survives for Rashba interactions a < 2\/ 3A (for chemical 
potential [i — 0), and has protected, helical edge states 
for geometries with open boundaries [31 [51 [UJ . We use 
t as the unit of energy (h = 1), take X/t = 0.2, and 
consider periodic lattices with L x L' unit cells. 

To investigate the effect of electron-electron repulsion, 
we consider the paradigmatic Hubbard interaction [42] 
and arrive at the Kane-Mele- Hubbard model 1431. 



KMH 



KM 



Hr 



H u = Wyp i c i -lf. (3) 



Hamiltonian ([3]) without Rashba coupling has been stud- 
ied intensely jl6 | 138 } [43l448| . In particular, its symmetries 
permit the application of exact quantum Monte Carlo 
methods without a sign problem [151 [Ml EH] • 

On a lattice with periodic boundaries, 7r fluxes can 
only be inserted in pairs, as illustrated for the minimal 
number of two fluxes in Fig. [Tj The flux pair is connected 
by a branch cut (or string), and every hopping process 
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spin fluxons (E = E n ) charge fluxons (E + , E~) 




FIG. 2. (Color online) In a quantum spin Hall insulator, a 
7r flux gives rise to four states (with charge q and spin S z ) 
localized near the flux, which lie inside the bulk energy gap 
between the valence and conduction bands (labeled "VB" and 
"CB" in the figure, respectively) |12l 131] . The states cor- 
respond to a Kramers doublet of spin fluxons |t) , |4-) with 
energy E^^, and a doublet of charge fluxons |+) , — ) with 
energies E + , E~ . 



crossing the cut acquires a phase e 17r = —1, as encoded by 
Tij in Eq. ([I]) . Different choices of the branch cut for fixed 
flux positions are related by a gauge transformation. 
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FIG. 3. (Color online) Spin susceptibility of the Kane-Mele 
model (X/t — 0.2) with two n fluxes at the maximal distance 
on an 18 x 18 lattice, for different Rashba couplings a. At 
temperatures ksT < O.lt, each tt flux contributes 2fc * T to the 
susceptibility, leading to \s ~ jr^f- Also shown is the spin 
gap energy scale 0.2A S for T — 0, a — 0. For a > 0, the 
chemical potential is adjusted to retain a half-filled band. 



IV. USING tt FLUXES TO PROBE 
CORRELATED TOPOLOGICAL STATES 



III. METHOD 

We have used the auxiliary-field quantum Monte Carlo 
method 03], which was previously applied to the Hub- 
bard model on the honeycomb lattice [15] , and the Kane- 
Mele-Hubbard model Q£l GHOSH]. The central idea of this 
stochastic method is to use a path integral representation 
of the interacting model ([3|. By means of a Hubbard- 
Stratonovich transformation, the Hubbard term is decou- 
pled, leading to a problem of noninteracting fermions in 
an external, space and imaginary-time dependent field. 
The sampling is over different configurations of these 
auxiliary fields in terms of local updates. For a given 
configuration of fields, Wick's theorem can be used to 
calculate arbitrary correlation functions from the single- 
particle Green function. We refer to a review [50], and 
previous work [15j [TBI HH] for technical details such as the 
calculation of energy gaps. 

Here we have used a projective formulation (with pro- 
jection parameter Ot = 40) to obtain ground state results 
starting from a trial wavefunction (the ground state of 
the U = case) 05], and a finite-temperature formula- 
tion to calculate thermodynamic properties. Both vari- 
ants rely on a Trotter discretisation of imaginary time 
(we used At = /3/L = 9/L = 0.1), but the associated 
systematic error is smaller than the statistical errors. At 
half filling, time-reversal invariance ensures that simula- 
tions can be carried out without a minus-sign problem 
even in the presence of tt fluxes. 



A. Thermodynamic signature of tt fluxes 

In the topological phase of the model ([I]), each tt flux 
gives rise to four fluxon states which are exponentially 
localized (due to the bulk energy gap A sp ) near the corre- 
sponding flux-threaded hexagons [12 GUJ , see Fig. [I] The 
states correspond to the spin fluxons |f) , ||) with energy 
E'^, forming a Kramers pair related by time reversal, 
and the charge fluxons |+) , |— ) (with energies E + ,E~) 
related by particle-hole symmetry. As we show in Fig. [2] 
the fluxon states lie inside the bulk band gap, and for 
noninteracting electrons = E + = E~ . 

The fluxons leave a characteristic signature in the 
static spin and charge susceptibilities 

Xs = (3 ((M z 2 ) - (M z ) 2 ) , Xc = (3 ((TV 2 ) - (AO 2 ) , 

(4) 

which are defined in terms of the operators of total spin 
in the z direction, M z = J2i^i aZ ^i' an d 01 the total 
charge, N — ^ c\c i ; the inverse temperature is given by 
P = jr^p- At low temperatures k^T <C A sp , we can re- 
strict the Hilbert space to {|f) , \i) , |+) , |— )}• If the spin 
fluxons are independent, and for a = 0, we expect a Curie 

law x s = Xc = 2feir P er n nux ' and hence Xs = Xc = k^r 
for two independent tt fluxes (see Appendix|A]). The pref- 
actor of the Curie law follows from the quantized spin 
Hall conductance in the absence of Rashba coupling [3] . 
Similarly, a Curie law was also predicted for topological 
excitations in polyacetylene [5T] . 

Figure [3] shows results for \s as a function of tem- 
perature for the Kane-Mele model with two tt fluxes 
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FIG. 4. (Color online) (a) Spin gap A s (q = 0) and single- 
particle gap A sp (q = K) in the thermodynamic limit as a 
function of the Hubbard repulsion U, at T = (A/t = 0.2, 
a = 0). (b) Scaling of S xx (L/2) using the critical exponents 
of the 3D XY model, z = 1, v = 0.6717(1) and /3 = 0.3486(1) 
|52| . The intersection gives the critical point U c /t = 5.70(3). 
The lattice size is L x L. Errorbars are smaller than the 
symbols. 
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FIG. 5. (Color online) Integrated dynamical spin structure 
factor Sn(i) at T — on a 9 x 9 lattice, (a) Localized spin 
fluxons created in the topological insulator phase at U/t — 4. 

(b) Absence of spin fluxons in the magnetic phase at U/t — 6. 

(c) Maximum of Sn(i), as a function of U/t. Here X/t — 0.2, 
a = 0, and Q/t = 0.2. 



located at the largest possible distance. At tempera- 
tures k^T « A s , Xs is dominated by bulk effects. For 
ksT < Q.lt, we observe the expected Curie law. The lat- 
ter is robust with respect to Rashba coupling, which is 
crucial for possible experimental realizations. 



B. Probing correlated topological insulators 

Figure [3] establishes the existence and thermodynamic 
signature of degenerate spin and charge fluxons in a quan- 
tum spin Hall insulator threaded by a pair of ir fluxes. 
We now consider the effect of electron- electron interac- 
tion in the framework of the Kane-Mele-Hubbard model 
([3]). For A > 0, the phase diagram of the latter includes 
a correlated quantum spin Hall insulating phase that is 
adiabatically connected to that of the Kane-Mele model 
(i.e., U = 0), and a Mott insulating phase with long- 
range antiferromagnetic order [131 EH]. Figure[4|a) shows 
the quantum phase transition between these two phases 
as a function of U/t at X/t = 0.2. At the transition, 
the spin gap A s — as obtained from finite-size scaling (see 
ref. 08] for details) — closes, corresponding to the conden- 
sation of magnetic excitons (47] 08] . The magnetic order 
is of the easy-plane type, and the transition has 3D XY 
universality corresponding to the ordering of local mo- 



ments 07] 08]- For U > U c , time- reversal symmetry is 
spontaneously broken, and the single-particle gap A sp 
remains open across the transition 05], see Fig. Qa). 

The location of the critical point can be estimated from 
the scaling behavior of the real-space spin-spin correla- 
tion function 

S xx (r) = (Sl(r)Sl(O)) (5) 

at the largest distance r = L/2. Here we consider 
correlations between spins on the A sublattice, but re- 
sults are the same for the B sublattice. The limit 
liniL^oo S xx (L/2) is identical to to 2 , with m being the 
magnetization per site. This critical value can be ob- 
tained by considering the 3D XY scaling behavior at the 
transition. Following ref.0SJ we plot L 2fs ^S xx (L/2) as a 
function of U for different system sizes using the critical 
exponents z = 1, v = 0.6717(1) and (3 = 0.3486(1) [55]. 
Figure |4|b) reveals the expected intersection of curves at 
the critical point, and gives U c /t = 5.70(3). 

The well-understood magnetic transition of the 
model |3]) provides a test case for the use of it fluxes to 
probe a correlated quantum spin Hall state, as well as to 
track the interaction-driven transition to a topologically 
trivial state. We solve the interacting model with two 
7T fluxes using exact quantum Monte Carlo simulations. 
Spin fluxons can be detected by calculating the lattice- 
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FIG. 6. (Color online) (a) Spin (xs) and charge (x c ) susceptibilities of the Kane-Mele-Hubbard model (X/t = 0.2, a = 0) at 
U/t — 4. We consider L x L lattices with one pair of tt fluxes placed at the maximal distance. At low temperatures, the spin 
susceptibility reveals a Curie law Xs = tt~t> whereas the charge susceptibility is suppressed by the charge gap. (b)-(c), Spin 
susceptibility as a function of temperature for different values of U/t (X/t = 0.2, a — 0). (a)-(c) show that with increasing 
U/t, the range of the interaction between spin fluxons increases, leading to deviations from the Curie law Xs = at low 
temperatures, (d) For U > U c = 5.70(3)t, Xs reflects the presence of long-range magnetic order in the bulk. Errorbars are 
smaller than the symbol size. The arrows indicate the energy scale associated with the spin gap. 



site resolved dynamical spin structure factor at T = 0, 
defined as 

S(i,u J )=nJ2\(n\cyti\V)\ 2 HEn-E -u J ) . (6) 

n 

Here AkmhI^) — E n \n), and |0) denotes the ground 
state. S(i,u>) corresponds to the spectrum of spin ex- 
citations at lattice site i. A real-space map of the spin 
fluxon states |j) , \\) is obtained by integrating S(i,oj) 
up to an energy scale fi/t = 0.2 well within the charge 

gap A c pa 2A sp , giving S n (i) = J Q dwS(i,u). For 
U/t = 4, corresponding to the quantum spin Hall phase 
[see Fig. Qa)], we see in Fig. [5^ a) very sharply defined 
spin fluxons localized at the two flux-threaded hexagons. 
The value of Sq(i) is about three orders of magnitude 
smaller at lattice sites which are further away from a 
flux, so that the spin fluxons can easily be detected nu- 
merically. In Fig. [5jb) , we show results for the magnetic 
insulating phase at U/t = 6. As expected for this topo- 
logically trivial state, no well-defined spin fluxons exists. 

The dependence of Sq(i) on U/t across the magnetic 
quantum phase transition is shown in Fig. [5^c). A clear 
signal is found deep in the topological insulator phase, 
whereas a strong drop is observed on approaching the 
critical point at U c /t — 5.70(3). Hence, the spin fluxon 
signal can be used in quantum Monte Carlo simulations 
to distinguish topological and nontopological phases. 

As for the noninteracting case (Fig. [3]), the spin flux- 
ons created by the tt fluxes give rise to a characteristic 
Curie law in the spin susceptibility. Figure |6]ja) shows 
quantum Monte Carlo results for the spin and charge 
susceptibilities defined in Eq. Q in the topological phase 
(U/t = 4). We again consider two tt fluxes at the max- 
imal separation. At low temperatures, k^/T -C A s , \s 
is well described by x s = or ttr P er 77 ^ ux - ^he 

additional factor of 2 compared to the case U = comes 
from the splitting of spin and charge states which only 



leaves the Kramers doublet {|f) , at low energies (see 
Appendix [A]) . The Curie law holds down to the lowest 
temperatures considered in Fig. [6]ja). Finally, the charge 
susceptibility is strongly suppressed at low temperatures, 
and reveals the absence of low-energy charge fluxons as 
a result of the Hubbard repulsion. 

C. Interaction between spin fluxons 

So far, we have exploited the thermodynamic and spec- 
tral signatures of independent spin fluxon excitations 
(i.e., free spins). On periodic lattices, spin fluxons can 
only be created in pairs, and it is therefore interesting to 
consider their mutual interaction. Such interactions will 
play a key role in Sec. |Vj where we consider quantum 
spin Hall insulators with multiple tt fluxes to create and 
simulate systems of interacting spins. 

Interaction effects due to a coupling between two spin 
fluxons in a lattice with one pair of tt fluxes become visi- 
ble for larger U/t, i.e., closer to the magnetic transition. 
Figures ^b) and (c) show a deviation from the Curie law 
below a temperature scale determined by the interaction 
between spin fluxons. In the model ([3]), this interaction 
is mediated by the exchange of collective spin excitations 
(magnetic excitons), which are the lowest lying excita- 
tions of the correlated topological insulator phase, and 
evolve into the gapless Goldstone mode of the magnetic 
state. Since magnetic order is of the easy-plane type, 
the dominant contribution of the resulting interaction is 
expected to have the general form 

S in t = -9 2 J2 ff dTdT ' (7) 
x [S+( T )D(r - r' ,t - t')S-,(t') + R.c] , 
where S^(t) are spin-flip operators acting on a spin 
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FIG. 7. (Color online) Spin susceptibility as a function of 
temperature for two tt fluxes arranged as shown in the in- 
set (X/t = 0.2, a = 0, 9 x 9 lattice). With increasing U/t, 
the strength of the interaction between spin fluxons increases, 
as revealed by the shift of the temperature below which de- 
viations from a t-^= Curie law occur. Statistical errors are 

k B T 

smaller than the symbol size. Inset: Sn(i) for U/t = 4, using 
the same color coding as in Fig. [5ja). 

fluxon at position r at time r, D(r, t) is the Fourier 
transform of the exciton propagator D(q,\Vl m ) (q: mo- 
mentum, Q m = 2nit/f3: bosonic Matsubara frequency), 
and g is a coupling constant. At long wavelengths, the 
dispersion relation of the collective spin mode can be 
written as uj(q) — \Jv 2 \q — Q\ 2 + Af, where v is the spin 
velocity, A s is the spin gap, and Q is the magnetic or- 
dering wavevector. The minimal exciton energy is given 
by u>{q = Q) = A s . Fourier transformation of the propa- 
gator (see Appendix [C]) gives in the limit of low energies 
and long wavelengths 

£(r,T)~exp(iQT)exp(-A s T)exp^-^i^ . (8) 

The first term determines the sign of the interaction. The 
decay at large imaginary time r is governed by the spin 
gap A s . The fast decay as a function of \r\ underlies 
the clear Curie law seen, e.g., in Fig. [6} The interaction 
range and strength can be tuned via the spin gap and 
hence [cf. Fi g ffl a)] by varying U/t. 

From Eq. (|8j) we expect the interaction range to in- 
crease with increasing U/t due to the decrease of A s , see 
Fig. [i|a). Indeed, Figs.^b) and (c) reveal an enhanced 
effect of the spin fluxon separation at low temperatures 
with increasing U/t. In particular, for U/t — 5.5 (close 
to the magnetic transition), Fig. [6^c) shows a Curie law 
corresponding to two free spin fluxons only for the largest 
system sizes (L = 18). As U — > U c , the interaction range 
diverges, and free spin fluxons can no longer exist. For 
U > U c , time-reversal invariance is broken and tt fluxes 
do not create spin fluxons. Instead, the spin susceptibil- 
ity in Fig. Rifd) is that of an antiferromagnet; the finite 
value of Xs/L 2 at T = reflects the density of spin wave 
excitations. 




k B T/t 

FIG. 8. (Color online) Spin susceptibility as a function of 
temperature for four tt fluxes arranged as shown in the inset 
(X/t = 0.2, a = 0, U/t = 4) on L x L lattices. The data 
reveal a Curie law at intermediate temperatures, and 

l^j; at low temperatures. Statistical errors are smaller than 
the symbol size. Inset: Sn(i) for L = 15, using the same color 
coding as in Fig. [5|a). 

To illustrate the dependence of the interaction strength 
on A s , we consider two fluxes at a fixed, small distance 
as illustrated in the inset of Fig. [7] We show the spin 
susceptibility for different values of U/t in Fig. [7j For 
U/t = 3, a Curie law Xs ~ may be inferred at tem- 
peratures k^T « O.li. Increasing U/t, the interaction 
between the spin fluxons becomes too large to observe 
free spin fluxons below the temperature range set by the 
bulk spin gap A s . The downturn of Xs occurs at higher 
and higher temperatures with increasing U/t, and reflects 
a tunable, correlation-induced energy scale for the inter- 
action between spin fluxons that is absent in Fig. [3j 

V. 7T-FLUX QUANTUM SPIN MODELS 

The possibility of inserting tt fluxes to create local- 
ized spin fluxons with a tunable interaction mediated by 
magnetic excitons provides a toolbox to engineer inter- 
acting spin models in correlated topological insulators. 
The computational effort for quantum Monte Carlo sim- 
ulations does not depend on the number of n fluxes, and 
the latter can be arranged in almost arbitrary geometries 
on the honeycomb lattice. 

A. Three-spin system 

As a first extension of the two-spin cases considered 
so far, we consider four spin fluxons emerging from two 
pairs of tt fluxes. The fluxes are arranged so that three 
spin fluxons form a ring, and the fourth spin fluxon is lo- 
cated at the largest distance from the center of the ring. 
For large enough lattices, the separated spin fluxon will 
not couple to the other three, and the physical problem is 
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FIG. 9. (Color online) Integrated local density of states 
An(i) (n = 0.2t, see text) at T = for the Kane-Mele model 
(X/t — 0.2, a = 0) with a periodic chain of n fluxes. We 
show a part of the 72 x 12 lattice used, and the size of the 
magnetic unit cell containing two n fluxes. The latter has 
width Orr = 4a, where a = 1 corresponds to the norm of the 
lattice vectors of the underlying honeycomb lattice. 



similar to experiments on coupled quantum dots [S3] or 
flux qubits |54j in the context of quantum computation. 
The three spin fluxons experience a transverse interac- 
tion of the form Q , and behave as an effective spin with 
S z = ±1/2 at low temperatures (see Appendi x [B| ). The 
spin susceptibility for U/t = 4 shown in Fig. IsJ reveals 
that, at low temperatures, the two independent spins 
indeed give rise to the expected Curie law Xs 
At higher temperatures k^T O.lt, we find Xt 
corresponding to four independent spin fluxons. In the 
regime where Xs = jr^f 3 the sign of the interaction be- 
tween the spin fluxons determines the ground state de- 
generacy of the three-spin cluster. A net ferromagnetic 
interaction results in a spin-1/2 doublet, whereas an an- 
tiferromagnetic coupling gives rise to a four-fold degen- 
erate, chiral ground state (see Appendix [Bj. In principle, 
the sign of the exchange coupling for the case of Fig.[8]can 
be determined from entropy measurements. Since Q = 
for the model Q, Eq. Q suggests that the interaction 
is ferromagnetic. 
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B. Simulation of one-dimensional fluxon chains 

Whereas the study of systems with a small number of 
spins is relevant for applications such as quantum com- 
puting, many questions in condensed matter physics are 
related to periodic spin lattices. In this section, we there- 
fore consider one-dimensional chains of tt fluxes in the 
honeycomb lattice with periodic boundary conditions. 

We begin with the noninteracting Kane-Mele model 
with a periodic flux chain. The fluxon excitations are vis- 
ible in Fig. [9] which shows the integrated local density of 
states Aq(i) — dujA(i, cj); the single-particle spectral 
function A(i, to) is defined as usual in terms of the single- 
particle Green function, A{i,(J) = -ImG(i,u). Whereas 
the fluxons are well localized in the direction normal to 
the chain, the overlap of neighboring fluxons in the chain 
gives rise to a tight-binding band inside the topological 
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FIG. 10. (Color online) Spectrum of eigenvalues of the Kane- 
Mele model with a periodic chain of 7r fluxes (cf. Fig. [9J . Here 
X/t = 0.2, a — 0, and the honeycomb lattice has dimensions 
72 x 12. Points correspond to eigenvalues, and lines to the 
band structure e(k) = ±2t sin(2fca) with t « 0.1264 and a = 1. 



band gap which can be seen in the spectrum shown in 
Fig. [10] The specific form of the band structure can be 
attributed to the fact that the smallest unit cell for the 
fluxon chain contains two flux-threaded hexagons (and is 
four hexagons wide) , see Fig. [9j Exploiting that the four 
possible fluxon states per hexagon, {|t) , \i) , |+) , |— )}, 
can formally be written in terms of the fermion Fock 
states , \i) , |0) , |t4)}j an d assuming nearest-neighbor 
hopping, a suitable Hamiltonian is given by 



H 



H.c. 



(9) 



where 4>, ip refer to the two flux-threaded hexagons in the 
unit cell, and i numbers the unit cells. The resulting band 
dispersion e(k) — ±2tsm(2ka) matches the low-energy 
bands in the spectrum (Fig. 10 ). The form of the effective 



low energy Hamiltonian and especially the gapless nature 
of the spectrum stems from the fact that the unit cell is a 
gauge choice; a translation by half a lattice vector, 
corresponds to a gauge transformation. This symmetry 
allows the intercell and intracell hopping integrals to dif- 
fer only by a phase e ie . Imposing time-reversal symmetry 
pins the phase factor to 9 = and 8 = tt, thus lead- 
ing to the dispersion relations ±2£cos [(k + 0/a 7r )a w /2]. 
The choice 8 = Tt produces the above mentioned disper- 
sion relation, and the choice 8 — merely corresponds 
to translating the reciprocal lattice by half a reciprocal 
lattice unit vector. 

In contrast to the helical edge states of a quantum spin 
Hall insulator, each of the two fluxon bands is spin de- 
generate. As a result, and because the system is half 
filled, we expect a Mott transition of charge fluxons for 



any nonzero electron-electron repulsion. Figure 11 shows 
the spin and charge susceptibilities of the Kane-Mele- 
Hubbard model on L x 12 lattices with L/2 tt fluxes and 
U/t = 4. The Hubbard U causes an exponential sup- 
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FIG. 11. (Color online) (a) Spin and (b) charge susceptibility 
of the Kane-Mele-Hubbard model (\/t = 0.2, a = 0) at U/t = 
4. We consider L x 12 lattices with L/2 7r fluxes arranged in 
a periodic one-dimensional chain. The inset in (b) shows the 
charge susceptibility as a function of inverse temperature on 
a logarithmic scale. The key in (a) applies to all panels. 



pression of the charge susceptibility at low temperatures 
[see Fig. |ll[ b) and inse t], w hereas low-energy spin fluxon 
excitations remain [Fig.[TTja)]. Hence, similar to the one- 
dimensional Hubbard model, the fluxon chain undergoes 
a Mott transition to a state with a nonzero charge gap 
but gapless spin excitations. 

In the Mott phase of the fluxon chain, the low-energy 
physics is expected to be described by spin fluctuations, 
and hence by an effective spin model with spins corre- 
sponding to Kramers doublets of localized spin fluxons. 
Because the interaction range depends exponentially on 
the spin gap, we expect nearest-neighbor interactions 
J xy , J zz between spin fluxons to dominate, except for 
the close vicinity of the magnetic transition. As argued 
before, the magnetic exciton is of predominantly easy- 
plane type, and we therefore expect anisotropic interac- 
tions, \J xy \ ^> \J ZZ \. The minimal model for the spin 
chain is the one-dimensional XXZ Hamiltonian, 

h = j» £ s ts*+i + J xy Y,( s ? s ^i + ■ (io) 

i i 

Using the ALPS 1.3 implementation [53] we can simu- 
late this model in the stochastic series expansion repre- 
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FIG. 12. (Color online) Spin susceptibility as a func- 
tion of temperature. Symbols correspond to quantum Monte 
Carlo results for the Kane-Mele-Hubbard model (X/t = 0.2, 
a — 0, U/t = 4) on L x 12 lattices with L/2 n fluxes ar- 
ranged in a periodic one-dimensional chain. Lines are quan- 
tum Monte Carlo results for the one-dimensional XXZ model 
with J zz /\J xy \ = -0.1 and L/2 lattice sites (spins). 



scntation to calculate the spin susceptibility as a func- 
tion of temperature. There is one free parameter, 
J xy /J zz , which is varied to obtain the best fit to the low- 
temperature susceptibility (at high temperatures, bulk 
states of the topological insulator begin to play a role) 
of the Kane-Mele-Hubbard model. For example, con- 
sidering six spins, a rather good match between the 
spin fluxon data and the XXZ model is obtained for 
J zz /\J xy \ = —0.1 (the sign of J xy is irrelevant), see 
Importantly, taking the same parameters, and 



Fig. 12 



simulating ten spins with both spin fluxons and the XXZ 
model, equally good agreement is found in Fig. |12| These 
results demonstrate that the spin fluxons form a one- 
dimensional spin system with well-defined interactions, 
and that a quantum spin Hall insulator with tt fluxes can 
indeed be used as a quantum simulator. 



VI. CONCLUSIONS 

In this work, we have presented quantum Monte Carlo 
results for a correlated quantum spin Hall insulator with 
topological defects in the form of tt fluxes. Such fluxes 
represent a universal probe for the topological index that 
can be used in the presence of electronic correlations, and 
does not rely on spin conservation or an adiabatic con- 
nection to a noninteracting topological insulator. Our 
results demonstrate that tt fluxes can be combined with 
exact numerical simulations, and lead to clear signatures 
of nontrivial topological properties in spectral and ther- 
modynamic properties. As a concrete example, we have 
studied the magnetic quantum phase transition of the 
Kane-Mele-Hubbard model at which time-reversal sym- 
metry is spontaneously broken. In principle, tt fluxes can 
also be used in connection with fractional quantum spin 
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Hall states. 

More generally, tt fluxes in correlated topological insu- 
lators allow one to construct and simulate quantum spin 
models, and hence lead to a novel class of quantum sim- 
ulators. This finding is not restricted to the Kane-Mele- 
Hubbard model considered here. In particular, mag- 
netism driven by electronic correlations — the origin of 
the interaction between spin fluxons — is a common phe- 
nomenon. The physics described here relies on the coexis- 
tence of magnetic correlations and time-reversal symme- 
try, and cannot be captured by static mean-field descrip- 
tions. The spin models share the topological protection 
of their host against, for example, disorder. In general, 
they are characterized by a dynamical, time-dependent 
interaction reminiscent of spin-boson problems. The de- 
tailed form and sign of the interaction, whose strength 
and range can be tuned via the electronic interactions, 
depends on the electronic Hamiltonian and the lattice 
geometry of the underlying topological insulator. Be- 
cause of the spin-orbit interaction, the spin symmetry is 
C/(l), and — similar to cold atom realizations of the quan- 
tum Ising model |56] — the spin-spin interaction is generi- 
cally anisotropic. We have provided explicit evidence for 
the feasibility of our idea in terms of simulations of two 
and four spins, as well as of one-dimensional spin chains. 
With additional Rashba terms, spin models with a dis- 
crete Z2 Ising symmetry result. Although spin fluxon 
states are still well defined [12 [3T] , it is a priori not clear 
which operators have to be measured in the numerical 
simulations. Finally, the concept of fluxons originating 
from tt fluxes carries over to three-dimensional topologi- 
cal insulators O E2] . 

An open question of central importance is whether 
the use of tt fluxes will enable us to study quantum 
spin systems which are currently not accessible to nu- 
merical methods, for example due to a sign problem in 
the presence of frustrated interactions. Whereas we have 
provided evidence for the possibility to simulate arrays 
and chains of quantum spins, and to tune the interac- 
tion strength and range, entropy measurements are re- 
quired to determine the sign of the interactions. How- 
ever, the latter are extremely demanding to carry out 
using discrete-time quantum Monte Carlo methods. A 
systematic effort to study spin fluxon chain and ladder 
geometries is currently in progress. 

Our idea may potentially also be used in experiments. 
A strongly correlated topological insulator on the honey- 
comb lattice may be realized with Na 2 Ir0 3 [34], or with 
molecular graphene [57]- It has been suggested that tt 
fluxes can be created in a quantum spin Hall insulator 
by means of an adjacent superconductor and a magnetic 
field [31]. This idea can be generalized to arrays of tt 
fluxes using Abrikosov lattices. Alternatively, tt fluxes 
may be realized using SQUIDs. A potential problem is 
that the diameter of the tt fluxes will not be of the order 
of the lattice constant. Other exciting recent proposals 
which are relevant for the realization of our idea include 
artificial semiconductor honeycomb structures [35] . cold 



atoms in optical lattices [12], and cold atoms on chips 
[58] . In solid state setups, tt fluxes can also be created 
by dislocations [35] [3U] or wedge disclinations [SS] . 
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Appendix A: Spin susceptibility for two tt fluxes 

A single tt flux in a topological insulator gives rise to 
four states, |f) , |L) , |+) , |— ). In the absence of correla- 
tions, these states are degenerate. At low temperatures, 
the spin susceptibility, defined in Eq. Q, can be calcu- 
lated using the Hilbert space formed by only these states. 
Defining an effective Hamiltonian = Ylih-^ IV 1 ) ("4>\ 
with tp e {+, -,f,l} and E+ = E~ = = E^ = E^, 
we obtain 

X S = P({M 2 Z )-{M Z ) 2 ) (Al) 
1 SyM^e-^-M _ 1 
k B T E^le-^M 2fc B T' 

For U >• k B T, the spin fluxons |t) , \i) are the only 
low-energy excitations, and \ s can be calculated by re- 
stricting ip to {t,|}- Since E^ — E^ — E^^ due to time- 
reversal symmetry, we get 

*- = i£r- (A2) 

For the case of two independent tt fluxes, the above 
results imply x s = j^t at J7 = 0, and x s = k~T ^ or 
U > 0. This agrees with the numerical results shown in 
Fig. [3] for [7 = 0, and in Figs. [§7] for U > 0. 

Our derivation is only valid in the absence of Rashba 
spin-orbit coupling a. However, the numerical results in 
Fig. [3] show that the low- temperature Curie law in x s is 
the same also for a/0. 

Appendix B: Spin susceptibility and ground state 
degeneracy for four n fluxes 

The results for the Kane-Mele-Hubbard model with 
four tt fluxes shown in Fig. |8] reveal a Curie law at 

low temperatures, and a -j^p Curie law at higher temper- 
atures. This finding can be understood as corresponding 
to either two or four noninteracting spins. The latter case 
corresponds to the spatially separate spin fluxon and an 
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effective spin- 1/2 Kramers doublet (formed by the three 
nearby spin fluxons) in the regime where x s ~ ^fr' an< ^ 
to four noninteracting spin fluxons in the regime where 

The cluster formed by the three nearby spin fluxons 
has the possible configurations 



\M Z \ = 



{|ttt),|IU}}, 



(Bl) 



\m z \ = \ : {\m) , \m) - \ut) , \m) , \m) , im» ; 



The eigenvalues are given by 

E(k) = 2Jcosk. 



(B4) 



For J < 0, the ground state has k = and energy 2J. For 
J > 0, the ground state is chiral, with k = ±27r/3 and 
energy —J. Taking into account the sector M z = —1/2, 
we find a total ground state degeneracy of two in the fer- 
romagnetic case J < 0, and fourm the antiferromagnetic 
case J > 0. 



where M z denotes the total spin in the z direction. 
Since the exciton-mediated interaction in the Kane-Mele- 
Hubbard model has the form given in Eq. (|7|) and hence 
promotes spin-flips, the ground state can be expected to 
have|AP| = 1/2. The above mentioned effective spin-1/2 
then corresponds to the two possible values M z = ±1/2. 

The degeneracy of the ground state depends on the 
sign of the interaction. Considering M z = 1/2, we have 
the allowed states \lfl), |t4t) and |ttl)- The spin- 
flip terms which connect these states are of the form 
J(S^~ +1 S~ + S~ +1 ) , with periodic boundary conditions. 
An equivalent representation is given by the Hamiltonian 



H=jj2(\j+m\ + \j){j+n) > 



(B2) 



which describes the hopping of a particle (the spin down) 
on a three-site ring, with |1) = |^tt) etc. The eigen- 
states are obtained by Fourier transformation, and have 
the form 



\k) 



1 - 



j) 



0,± 



27T 



(B3) 



Appendix C: Fourier transform of the exciton 
propagator 



The exciton propagator in Eq. ^ takes the form 



D(q,r) = 



-™(q) 



e Pu(q) — 1 g-/30)(g) _ I 

with the exciton energy 



(CI) 



a;(g) = v /« 2 |g-Q| 2 + A2. (C2) 
In the low-temperature limit /3 — > oo, the propagator 
becomes D(q,r) = e~™ {q) . Setting q' = q - Q, the 
Fourier transform is given by 



D{r,T) = e iQ r J d 2 q'e iq ' r e-^ q ' +Q ^ T 
Assuming A s 3> v, we can expand to obtain 



v 2 \q'\ 2 
2A? 



(C3) 



(C4) 



Taking the continuum limit, the Fourier transform in- 
volves the product of two Gaussian integrals, and the 
result is given by Eq. ((8|). 
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